In [1]:
import visual_behavior_glm.src.GLM_params as glm_params
import visual_behavior_glm.src.GLM_analysis_tools as gat
from visual_behavior_glm.src.glm import GLM
import matplotlib.pyplot as plt
import visual_behavior.data_access.loading as loading
import visual_behavior.database as db
import plotly.express as px
import pandas as pd
import numpy as np
import os
from sklearn.cluster import KMeans
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
import seaborn as sns
In [2]:
%matplotlib inline
get results for glm version 4¶
In [3]:
results_all = gat.retrieve_results(results_type='full')
In [4]:
results_all['glm_version'].unique()
#use v4
Out[4]:
In [5]:
rs = gat.retrieve_results(search_dict = {'glm_version': '4_L2_optimize_by_cell'}, results_type='summary')
len(rs)
Out[5]:
make a new column called 'identifier' that is {OPHYS_EXPERIMENTID}{CELL_SPECIMEN_ID}. This will give each a unique ID¶
In [6]:
rs['identifier'] = rs['ophys_experiment_id'].astype(str) + '_' + rs['cell_specimen_id'].astype(str)
In [7]:
rs
Out[7]:
pivot the table to get it into a form for clustering¶
get fraction change in explained variance after dropout of each regressor¶
In [8]:
model_output_type = 'fraction_change_from_full'
rsp = rs.pivot(index='identifier',columns='dropout',values=model_output_type).reset_index()
rsp
Out[8]:
get variance explained for full model¶
In [9]:
tmp = 'variance_explained'
ve = rs.pivot(index='identifier',columns='dropout',values=tmp).reset_index()
ve
Out[9]:
In [11]:
tmp = ve.rename(columns={'Full':'varience_explained_full_model'})
rsp = rsp.merge(tmp[['identifier','varience_explained_full_model']], on=['identifier'])
filter out poorly fit cells¶
In [10]:
cells_to_include = ve[ve['Full']>0.01].identifier.values
order = np.argsort(ve[ve.identifier.isin(cells_to_include)==True]['Full'])
cell_order = cells_to_include[order]
len(cells_to_include)
Out[10]:
In [12]:
rsp = rsp[rsp.identifier.isin(cells_to_include)==True]
len(rsp)
Out[12]:
merge in the identifier, cre_line, session_type and equipment_name columns¶
In [13]:
rspm = rsp.merge(rs[['identifier','cre_line','session_type','equipment_name']].drop_duplicates(),left_on='identifier',right_on='identifier',how='inner')
rspm
Out[13]:
add a column for session ID (1, 3, 4, or 6)¶
In [14]:
def map_session_types(session_type):
session_id = session_type[6:7]
return session_id
rspm['session_id'] = rspm['session_type'].map(lambda st:map_session_types(st))
rspm['session_id'].unique()
Out[14]:
save feature matrix to disk¶
In [15]:
# save = True
# if save:
# rspm.to_csv('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/fraction_change_var_explained_v_4_L2_fixed_lambda=1_filtered_2020.08.11.csv', index=False)
get the column names to use for clustering¶
In [16]:
cols_for_clustering = [col for col in rspm.columns if col not in ['identifier','cre_line','session_type','equipment_name', 'session_id']]
cols_for_clustering = [col for col in cols_for_clustering if col not in ['image0','image1','image2','image3','image4','image5','image6','image7',
'visual']]
cols_for_clustering
Out[16]:
In [17]:
cols_for_clustering = [
'omissions',
'all-images',
'image_expectation',
'change',
'hits',
'misses',
'correct_rejects',
'false_alarms',
'post_lick_bouts',
'post_licks',
'pre_lick_bouts',
'pre_licks',
'rewards',
'pupil',
'running',
'time',
'model_bias',
'model_omissions1',
'model_task0',
'model_timing1D',
]
plot feature matrix¶
In [18]:
feature_matrix = rspm[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, cmap='RdBu_r', ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[18]:
In [19]:
# sorted by overall variance explaned
feature_matrix = rspm.sort_values('varience_explained_full_model').reset_index()[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, cmap='RdBu_r', ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[19]:
In [20]:
# sorted by omission dropout
feature_matrix = rspm.sort_values('omissions').reset_index()[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, cmap='RdBu_r', ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[20]:
In [ ]:
do k-means on the full feature matrix¶
In [21]:
# create kmeans object
kmeans = KMeans(n_clusters=4)
In [22]:
# fit kmeans object to data
kmeans.fit(rspm[cols_for_clustering])
Out[22]:
In [23]:
rspm['k-means_cluster_id'] = kmeans.fit_predict(rspm[cols_for_clustering])
how many clusters do we get?¶
In [24]:
rspm['k-means_cluster_id'].value_counts()
Out[24]:
make a tsne plot¶
In [25]:
tsne = TSNE(n_components=3, verbose=1, perplexity=40, n_iter=300)
In [26]:
tsne_results = tsne.fit_transform(rspm[cols_for_clustering])
assign tsne reults to new columns¶
In [27]:
rspm['tsne-2d-one'] = tsne_results[:,0]
rspm['tsne-2d-two'] = tsne_results[:,1]
In [28]:
rspm['session_type'].unique()
Out[28]:
make a histogram of the omission dropout values¶
In [29]:
fig,ax=plt.subplots()
ax.hist(rspm['omissions'],bins=np.arange(-1,1,0.01))
ax.axvline(0,color='black')
ax.set_xlabel('omissions dropout fractional change in var explained')
Out[29]:
make a new column that identifies cells as 'modulated' by running or omissions (the threshold is up for debate!)¶
In [30]:
def is_modulated(v, threshold=0):
return v < threshold
rspm['omission_modulated'] = rspm['omissions'].map(lambda v:is_modulated(v, threshold=0))
rspm['running_modulated'] = rspm['running'].map(lambda v:is_modulated(v, threshold=0))
plot the tsne results, coloring dots in different ways¶
It looks like they cluster almost perfectly by omission modulation!
In [31]:
fig,ax=plt.subplots(2,2,figsize=(15,12),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="cre_line",
palette=sns.color_palette("hls", 3),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[0,0]
)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_type",
palette=sns.color_palette("hls", 8),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[0,1]
)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="k-means_cluster_id",
palette=sns.color_palette("hls", 4),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[1,0]
)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue='omission_modulated',
palette=sns.color_palette("hls", 2),
data=rspm,
legend="full",
alpha=0.3,
ax=ax[1,1]
)
fig.tight_layout()
for ii,hue in enumerate(['cre_line','session_type','k-means label', 'omission modulated']):
ax.flatten()[ii].set_title('color by {}'.format(hue))
fig.tight_layout()
plot tsne again, grouped by cre-line¶
In [32]:
fig,ax=plt.subplots(2,2,figsize=(20,12),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm, #.query('cre_line=="Vip-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[0,0]
)
ax[0,0].set_title('all cre lines')
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Vip-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[0,1]
)
ax[0,1].set_title('VIP only')
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Slc17a7-IRES2-Cre"'),
legend='full',
alpha=0.3,
ax=ax[1,0]
)
ax[1,0].set_title('Slc17a7 only')
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Sst-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[1,1]
)
ax[1,1].set_title('Sst only')
Out[32]:
plot tsne again, colored continuously by the omissions dropout¶
In [33]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue='omissions',
palette='magma',
data=rspm.query('omissions>-1 and omissions<1'),
legend=False,
alpha=0.3,
ax=ax
)
Out[33]:
do kmeans clustring on the tsne space¶
In [34]:
kmeans = KMeans(n_clusters=2)
kmeans.fit(rspm[cols_for_clustering])
tsne_cols = ['tsne-2d-one','tsne-2d-two']
rspm['k-means_on_tSNE'] = kmeans.fit_predict(rspm[tsne_cols])
rspm['k-means_on_tSNE'].value_counts()
Out[34]:
In [35]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x="tsne-2d-one",
y="tsne-2d-two",
hue='k-means_on_tSNE',
palette=sns.color_palette("hls", 2),
data=rspm,
legend='full',
alpha=0.3,
ax=ax,
)
Out[35]:
UMAP¶
In [105]:
import umap
min_dist = 0
spread=5
n_neighbors = 10
metric = 'euclidean'
reducer = umap.UMAP(n_components=2, n_neighbors=n_neighbors, min_dist=min_dist, metric=metric, spread=spread)
mapper = reducer.fit(rspm[cols_for_clustering])
embedding = reducer.fit_transform(rspm[cols_for_clustering])
embedding.shape
Out[105]:
In [106]:
mapper.embedding_.shape
Out[106]:
In [107]:
# save_dir = r'\\allen\programs\braintv\workgroups\nc-ophys\visual_behavior\clustering\UMAP_output'
In [108]:
# filepath = os.path.join(save_dir, '200810_UMAP_10D_GLM_fraction_change_in_explained_variance')
# np.save(file=filepath, arr=embedding_3d)
In [109]:
# filepath = os.path.join(save_dir, '200810_UMAP_3D_GLM_fraction_change_in_explained_variance')
# np.save(file=filepath, arr=embedding_3d)
In [110]:
# save_dir = r'\\allen\programs\braintv\workgroups\nc-ophys\visual_behavior\clustering\UMAP_output'
# filepath = os.path.join(save_dir, '200810_UMAP_2D_GLM_fraction_change_in_explained_variance.npy')
# data = np.load(filepath)
assign umap results to new columns¶
In [111]:
rspm['umap_embedding_0'] = embedding[:, 0]
rspm['umap_embedding_1'] = embedding[:, 1]
# rspm['umap_embedding_2'] = embedding[:, 2]
In [112]:
umap_cols = ['umap_embedding_0', 'umap_embedding_1']#, 'umap_embedding_2']
do kmeans on the 2d embedding¶
only two clusters since the 2d umap gives two clear clusters k-means seems to identify them perfectly!
In [113]:
kmeans = KMeans(n_clusters=3)
umap_cols = ['umap_embedding_0','umap_embedding_1']
rspm['k-means_on_umap'] = kmeans.fit_predict(rspm[umap_cols])
rspm['k-means_on_umap'].value_counts()
Out[113]:
plot the results in 2D¶
k-means clearly didn't do a good job of identifying two of the three clusters
In [114]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
sns.scatterplot(
x='umap_embedding_0',
y='umap_embedding_1',
hue='k-means_on_umap',
palette=sns.color_palette("hls", 3),
data=rspm,
legend='full',
alpha=0.3,
ax=ax,
)
Out[114]:
In [115]:
import seaborn as sns
fig,ax=plt.subplots(2,2,figsize=(15,12),sharey=True,sharex=True)
sns.scatterplot(
x='umap_embedding_0',
y='umap_embedding_1',
hue="cre_line",
palette=sns.color_palette("hls", 3),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[0,0]
)
sns.scatterplot(
x='umap_embedding_0',
y='umap_embedding_1',
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[0,1]
)
sns.scatterplot(
x='umap_embedding_0',
y='umap_embedding_1',
hue='k-means_on_umap',
palette=sns.color_palette("hls", len(rspm['k-means_on_umap'].unique())),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[1,0]
)
sns.scatterplot(
x='umap_embedding_0',
y='umap_embedding_1',
hue='omission_modulated',
palette=sns.color_palette("hls", 2),
data=rspm,
legend="full",
alpha=0.1,
ax=ax[1,1]
)
fig.tight_layout()
for ii,hue in enumerate(['cre_line','session_type','k-means label', 'omission modulated']):
ax.flatten()[ii].set_title('color by {}'.format(hue))
fig.tight_layout()
In [116]:
fig,ax=plt.subplots(2,2,figsize=(20,12),sharey=True,sharex=True)
sns.scatterplot(
x="umap_embedding_0",
y="umap_embedding_1",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm, #.query('cre_line=="Vip-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[0,0]
)
ax[0,0].set_title('all cre lines')
sns.scatterplot(
x="umap_embedding_0",
y="umap_embedding_1",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Vip-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[0,1]
)
ax[0,1].set_title('VIP only')
sns.scatterplot(
x="umap_embedding_0",
y="umap_embedding_1",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Slc17a7-IRES2-Cre"'),
legend='full',
alpha=0.3,
ax=ax[1,0]
)
ax[1,0].set_title('Slc17a7 only')
sns.scatterplot(
x="umap_embedding_0",
y="umap_embedding_1",
hue="session_id",
palette=sns.color_palette("hls", 4),
data=rspm.query('cre_line=="Sst-IRES-Cre"'),
legend='full',
alpha=0.3,
ax=ax[1,1]
)
ax[1,1].set_title('Sst only')
Out[116]:
now make a 3d plotly plot of the same¶
In [117]:
import plotly.express as px
fig = px.scatter_3d(
rspm,
x='umap_embedding_0',
y='umap_embedding_1',
z='umap_embedding_2',
color='k-means_on_umap'
)
fig.show()
plot umap with omissions modulation as continuous color¶
In [118]:
fig,ax=plt.subplots(1,1,figsize=(10,6),sharey=True,sharex=True)
ax = sns.scatterplot(
x='umap_embedding_0',
y='umap_embedding_1',
hue='omissions',
palette='magma',
data=rspm.query('omissions>-1 and omissions<0.25'),
legend=False,
alpha=0.3,
ax=ax,
)
# plt.colorbar(cax=ax, title='omission_modulation')
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: